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Abstract 

We propose in this paper a random intercept Poisson model in which the random effect distribution 
is assumed to fohow a generahzed log-gamma (GLG) distribution. We derive the first two moments 
for the marginal distribution as well as the intraclass correlation. Even though numerical integration 
methods are in general required for deriving the marginal models, we obtain the multivariate negative 
binomial model for a particular parameter setting of the hierarchical model. An iterative process 
is derived for obtaining the maximum likelihood estimates for the parameters in the multivariate 
negative binomial model. Residual analysis are proposed and two applications with real data are 
given for illustration. 

Key words : Count data; Generalized log-gamma distribution; Multivariate negative binomial distri- 
bution; Overdispersion; Random-effect models. 

1 Introduction 

The effects of the misspecification of the random effect distribution in generalized linear mixed models 
(GLMMs) have been investigated by some authors recently. For instance, Litiere et al. (2008) verified 
by Monte Carlo studies that the misspecification of the random effect distribution of the response 
variable in random intercept logistic models may lead to severe bias in the random effect component 
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prediction, which in many problems may be the main parameter of interest. These same authors have 
proposed a family of tests to detect the misspecification of the random effect distribution in GLMMs 
(Alonso et al., 2008). In addition, Lcc and Ncldcr (1996, 2001) have suggested a flexibilization of the 
random effect distribution in GLMMs, but under a hierarchical framework. Although any combination 
between the conditional response and the random effect distributions may be considered in the Lee 
and Nelder's proposal, the majority of the applications have been done for conjugate distributions. In 
particular, under the marginal framework, Molenberghs et al. (2007) have presented a combination 
between gamma and normal random effects in Poisson mixed models and more recently Zhang et al. 
(2008) assumed a log-gamma distribution for the random effects in linear mixed models. 

The aim of this paper is to present an alternative distribution for the random effect in random 
intercept Poisson models, which is characterized by assuming a generalized log-gamma distribution for 
the random effect component. This distribution introduced by Prentice (1974) (see also Lawless, 1980) 
has as particular cases the normal and extreme value distributions and it may assume skew forms to 
the right and to the left. In addition, generalized log-gamma models have been widely applied in the 
areas of survival analysis and reliability. For instance, DiCiccio (1987) derived approximate inferences 
for the quantiles and scale parameters whereas Young and Bakir (1987) obtained the bias of order n~^, 
where n is the sample size, for the parameter estimates in generalized log-gamma regression models for 
uncensored samples. Young and Bakir (1987) also presented the expectation of various log-likelihood 
derivatives in closed-form expressions. Ahn (1996) proposed a regression tree method to classify 
the heterogeneous subsets of the data into different generalized log-gamma regression models with 
the shape parameter being estimated separately in each formed stratum under independent random 
censoring. Ortega et al. (2003) derived the appropriate matrices for assessing local influence on the 
parameter estimates under different perturbation schemes and Chien-Tai et al. (2004) presented a 
conditional method of inference to derive confidence intervals for the location as well as quantiles and 
reliability functions under progressively type-II censoring and by assuming the shape parameter known 
in generalized log-gamma regression models with censored data. More recently, Cox et al. (2007) 
presented a taxonomy of the hazard function of generalized gamma distribution with application to 
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study of survival after diagnosis of clinical AIDS during different phases of HIV therapy and Ortega et 
al. (2009) introduced the generalized log-gamma regression models with cure fraction giving emphasis 
to assessment of local influence. 

The paper is organized as follows. In Section 2 wc present a brief review of the generalized 
log-gamma distribution. The random intercept Poisson generalized log-gamma model is proposed in 
Section 3, as well as a discussion on the parameter and random effect estimation. The derivation of the 
first two moments for the marginal distribution and of the intraclass correlation are given in Section 4. 
For a particular parameter setting of the hierarchical model we derive, in Section 5, the multivariate 
negative binomial model (Johnson et al., 1997) as a marginal model. An iterative process for the 
parameter estimation as well as goodness-of-fit procedures and residual analysis are also presented in 
Section 5. In Section 6 the epilepsy data set (Diggle et al., 2002) is fitted with the random intercept 
Poisson-normal and random intercept Poisson-GLG models and compared under the AIC criterion. 
Another application that has been analyzed by Poisson models (Lange et al., 1994) is reanalyzed with 
the negative binomial and multivariate negative binomial models. 

2 Generalized log-gamma distribution 

Let y be a random variable following a generalized log-gamma distribution. The probability density 
function (pdf) of y (see, for instance. Lawless, 2002) is given by 



where y E SI; /i E M, a > and A G iR are, respectively, the position, scale and shape parameters 



The extreme value distribution is a particular case of (1), when A = 1. For A < the pdf of y is 
skew to the right and for A > it is skew to the left. Figure 1 presents the behavior of the pdf of 
y GLG(0, 1, A) for some values of A. 




(1) 



and c(A) 



■^{X '^)^ with r(-) being the gamma function. We will denote y ^ GLG(//, a. A). 
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One has for A 7^ the following moments: 

E{y) = fi + al^ I and Var(y) = \ 

where 'ip{-) and denote, respectively, the digamma and trigamma functions. For A = one has 
the normal case for which E(y) = ^ and Var(y) = a^. 
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Figure 1: Graphs of the generalized log-gamma distribution for some values of A and by assuming 
jjL = Q and £7 = 1. 

3 The random intercept Poisson-GLG model 

Let yij denote the jth outcome measured for the ith cluster (subject), i = 1, . . . , n and j = 1, . . . , mj. 
We will assume the following random intercept Poisson-GLG model: 

(i) yij\bi Poisson(nij) 

(ii) Uij = exp(xT./3 + 6j) and 

(iii) bi - GLG(0,c7,A), 



X = 2 

X^-Z 
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where Xy = {xiji, . . . ,Xijp)~^ contains values of explanatory variables and /3 = . . . ,Pp)~^ is the 
parameter vector of the systematic component. The model (i)-(iii) will be named random intercept 

Poisson-GLG model. When A = one has the random intercept Poisson-normal model (see, for 
instance, Brcslow and Clayton, 1993). Let fY\biyij\bii P) and ff,(bi;a,X) be the probability function 
of yij\bi and the pdf of bi, respectively. Then, the marginal probability function of y — (yij • • • > Yn) > 
where = (ya, ■ ■ ■ , VimiV , is given by 

n 

fY{y;l3,a,X) = l[fYiyi;f3,a,X) 

i=l 

with 

f + OO ( "^i 1 

fY{yi;P,a,\) = l]JfY\biyij\h,P)}Mbi;a,X)dbi, (2) 

which in general does not have a closed-form. Then, the log-likelihood function for the marginal 
model, using (2), takes the form 

n I.+00 ( ] 

Lie) = J2^og m fYibiyij\bi,f3) \ hibi;a,X)dbi, (3) 

i=l J-oo J 

where 6 = (/3^,A,(t)^. Expression (3) should be approximated by numerical integration methods, 
such as Laplace approximation or Gauss-Hermite quadrature. To predict the random effects we can 
use the empirical Bayes method (see, for instance, McCulloch and Searle, 2001) given by 

b - E[b \y ] - ^^^^^ 
J^^ fY\b{yi\bi,fB)fb{bi;a,X)dbi 

The NLMIXED procedure available in SAS has been required to obtain the parameter estimate in the 
GLMM class. Through this procedure it is possible to compute the integral in (3) and to perform the 
maximization of the log likelihood in (3) as well as to obtain the random effect prediction. 

4 Derivation of moments 

We derive in this section the moments E(yij) and Var(j/jj) as well as Cov(yij,yiji), for j / j' , and the 
following results will be used: 



5 



(a) Eivij) = E{E{yij\bi)} = //iiE(e^'), 

(b) Var(j/i,) = Mar{E{yij\bi)} + E{V&v{yij\bi)} = 4Var(e*'0 + /^i,E(e^O and 

(c) Cov{yij,yif) = Cov{E{yij\bi),Eiyiji\bi)} + E{Cov{yij,yijf\bi)} = 

= CoY{iiije^i,iiijie^i) + = iXijiiij' Var(e*0> for 3 + /> 

where /Xy = exp(x^/3). From (a)-(c) above one has that 

Var(j/tj) _ Var(e^') 

and since ix^ > 0, Var(e''*) > and E{e'^^) > it follows that Var(yjj) > E{yij), that is, the model 
(i)-(iii) is overdispersed. In Table 1 one has the expressions of E{e^^) and E(e^''') for some ranges of 
A, where 

poo „ poo „ 

h{X,a)= / t^-'^^'^+^y^e-'dt and h{\,(T) = / (^Aa+i)-!^-*^^^ 
Jo Jo 

which should be solved by numerical integration methods. 

Table 1 



First two moments for the random variable e**' according to the values of the 

shape parameter A. 



Shape parameter 


Eie"*) 


Eie^"^) 


A = 






A > 
A < 


0^r{A-(Aa + l)} 


r(x-2) l2(A,o-) 



5 The multivariate negative binomial model 

Consider now the following random intercept Poisson-GLG model: 

(i) Uijlh ™' Poisson(tiij), 

(ii) Uij = exp(x^/3 + bi) and 

(iii) bi - GLG(0,A,A), 
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that is, the same model given in Section 3 with a = X (A > 0). 
Denoting = the joint distribution of (yj, bi) is given by 

{mi 1 
n fY\b{yij\bi, ^) I fb{bi;(l)) 



where = Yl^=i Vij ^^id = X^^^ fJ-ij. Consider the variable transformation ti = exp(6j)(/Xi+ + ^). 
One has that ^ = j: and the joint distribution of {yi,ti) assumes the form 



r(</>)(nf=iy.,!)(</' + m)(<^+^^+) 

Thus, the marginal probability function of yields 

I mi Vij^ 
j=l H'ij , 

^(<^)(^^il^/^,■!)(</> + /^^+)(*+^^+) 

and since T{(f) + = /q°° e~^Hf~^'^'^~^dti, the marginal probability function of y^ reduces to 

that is the multivariate negative binomial distribution (see, for instance, Johnson ct al., 1997) of means 

2 

E{yij) = fiij, variances Var(y.y) = fiij + for j = 1, . . . ,mj, and covariances Co\ {yij,yiji) = , 
for j / /. The intraclass correlation between y^j and yjj/, for j / /, can be expressed as 

CoTT{yij,yif) = " 

These correlations are always positive and for large values of (f) the multivariate negative binomial 
counts y'^jS behave approximately as independent Poisson observations with respective means iJ,[jS. 

Therefore, we derive the multivariate negative binomial distribution from an alternative way, by 
assuming a particular log-gamma distribution for the random effect in a random intercept Poisson 
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model. By a similar calculation one may show that the marginal distribution of yij is a negative 
binomial distribution of mean ^ij, variance /i^ + and dispersion parameter > (see, for instance, 
McCullagh and Ncldcr, 1989). We will denote ~ MNB(/Xj, 0), where y^ = {yn , ■ ■ ■ , yirrH)~^ , fJ-i = 
{jjLii, . . . , fiirni)^ and > 0. The multivariate negative binomial model is defined by assuming that 
lo g/^ij — Then, the log-likelihood function for the multivariate negative binomial model yields 

i=i ^ ^\^) J j=i j=i i=i j=i 

n rrii I xJ.d 



+EE^^>g | ^^^^ h (5) 

i=ij=i \(p + l2j=ie J 

where 6> = (/S"^, 0)"^. The score function and the Fisher information matrix may be obtained for 9 
and an iterative process can be performed to get the maximum likelihood estimates. 

5.1 Score function 

The score function is obtained from (5) by derivating the log- likelihood function L(0) with respect to 
/3 and 0, respectively. We obtain 

n 

= E^nyi-«i/*i) and 

i=l 

by using the result r(^ -|- yi+)/T{(l)) = (f){^ + l){<p + 2) ...{(/) + yi+ — 1) (see, for instance. Lawless, 
1987) we find 

n f(2/<+-l) 

yi+ i„„n , ^-1.. ^ , /^iH 



u<^ = E E (j + 0)-' - - Ml + r Vi-h) + 



j=l j=o VV + W+) (^ + W+) 

where Oj = , Xj is an mj x p matrix of rows x^, for j = 1, . . . , mj, so that J2jtQ~^^ is zero 

when — 1 < 0. 
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5.2 Fisher information matrix 

The Fisher information matrix for /3 is obtained such that 



K 



/3/3 





dV 


[ dp 








[s[ 






n 




E 

i=l 


x7{ 



/3 



yf_ + yi+)_-S~^ w T yi+) \ ^ \ ' 

i=l 

in that, D(/ij) = diag {^Uji, . . . , /t/jmi} • The calculation of the Fisher information for (j) follows the 
similar steps of the negative binomial model (see, for instance. Lawless, 1987). We find 



(0 + yi+) 



mi 



EE 



=1 

oo 



E (j+^y 

j=0 



yi+ 



E|E(i + <^)"'P(^^+>J') 

i=0 [j=0 



A-1 



+ 



In addition, it may be showed the orthogonahty between /3 and 4>, as in the univariate case. Thus, 
the Fisher information matrix for takes the block-diagonal form Kqo = diag{K^^, K^^}. 

5.3 Iterative process 

Similarly to the univariate case we can perform a scoring Fisher and a Newton-Raphson iterative 
processes for obtaining the maximum likelihood estimates /3 and (p, respectively, which are given by 



(0) 



and 



V^, for r = 0,l,2,.... 



(7) 
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in that Wj = D(/ij) — {cj) + fiij^)^"^ ji^iij and L^,^ = dV^/dcj). We can start the iterative process defined 
in ([6]) and ([7]) by using, for instance, the maximum hkehhood estimates from the univariate case in 
which the rrii observations for each group are assumed independent. For large sample (n large) we 
expected that the maximum likelihood estimators follow, under suitable regularity conditions, normal 
distributions. That is, for n large /9 ~ Np(/3,K^^) and (f> ~ N((;^), K^J). In addition, due to the 
orthogonality between (3 and (j) one has the asymptotic independence between (3 and (f). 

5.4 Residual analysis 

We found that the estimates of the multivariate negative binomial saturated model log-likelihood 
function is fi^j = yij, and therefore, the MNB deviance function has the following expression, 



Z?(y,A„</') = E2 



i=l 



^^^^ \ ir> r + ' 



(8) 



Similarly to the univariate case (see, for instance, Svetliza and Paula, 2003) we can define the de- 
viance as a measure for multivariate negative binomial models from ([S]). So, after some algebraic 
manipulation and by assuming that (j) is fixed and yij > 0, Vzj, we may express the deviance as 
D(y; (i) = Ya=i Ej^i d'^(.yij,f^ij,(p), where 

The quantity (p that appears in the deviance expression may be replaced, for instance, by a consistent 
estimate of <j), such the maximum likelihood estimate (p. 

As a residual proposal we can work, for instance, with the deviance component residual similarly to 
the univariate case. Svetliza and Paula (2003) performed various simulation studies with indication of 
a very good agreement between the empirical distribution of the deviance component residual and the 
normal distribution, even for (j) small. In the multivariate case we will adopt the following expression 
for the deviance component residual: 



d*iyij; fj-ij,(l)) 
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where 

in that {yij , fuj , (/>) is defined in ([9]) , hijj is the jth principal diagonal element of the matrix Hj = 
'Wy^X.i{'X.J'WiX.i)~^X.J'Wy^ and the sign is the same of [yij — jiij). 

A suggestion in order to assess departures from the postulated distributions for the responses 
y'ljS, as well as the presence of outlying observations, is to perform the normal probability plot for 
d{yij , jlij , (j)) with generate envelope (see, for instance, Svetliza and Paula, 2003). Another possibility, 
suggested by Waller and Zelterman (1997), is the Pearson residual whose expression is given by 

— / 

Again, here it is recommended to perform the normal probability plot with generated envelope to 
detect possible departures from the error assumptions as we as outlying observations. 

6 Applications 
6.1 Epilepsy data 

Diggle et al. (2002) described an experiment in which 59 epileptic patients were randomly assigned to 
one of two treatment groups: treatment (progabide drug) and placebo groups. The number of seizures 
experimented by each patient during the baseline period (eight-week) and the four consecutive periods 
(two-week) was recorded. The main objective of this application is to analyze the drug effect and 
compare its effect with the placebo group effect. Overdispersion evidences were observed in the data 
set and a generalized estimating equation model was applied to fit the data. In order to illustrate the 
potentiality of the random intercept Poisson-GLG model we modify the data set, patient #49 was 
dropped and some values of the patient #18 were modified in order to make it an outlying observation. 
Then, similarly to Diggle et al. (2002), we assume the following random intercept Poisson-normal 
model: 

(i) y-ikjlk '~ P(iijfcj), 



11 




8 10 12 14 16 

Time 



Figure 2: Profile of the seizures of the placebo and treatments groups. 

(ii) Uiio = cxp(a + k) + log (to), 

Uiij = cxp{a + /3 + bi) + log{tj), i = 1, . . . , 28, 
Ui20 = exp(Q + bi) + log (to), 

Ui2j = exp(a + P + S + bi) + log(tj) i = 29, . . . , 59, and 

(iii) 6,'^^N(0,a2), 

in that, pikj denotes the number of seizures experienced by the ith patient in the kth. group and _7th 
period, where z = 1, . . . , 59, A; = 1, 2 (placebo(28) and treatment (31)) and j = 1,2,3,4. In addition, 
tj denotes the week number of the jth period {to = 8 and ti = t2 = ts = t4 = 2), /3 is the parameter 
which represents the treatment effect and S is the parameter referents to the treatment group effect 
in relation to the placebo group. 

However, from Figure 2, that describes the log(counting) of seizures in the baseline period, we 
notice a skew form to right suggesting a skew distribution for the random intercept. Thus, a random 
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intercept Poisson-CLG model (with A < 0) is also assumed to fit this data set. Indeed, we replace 
in the model (i)-(iii) the assumption bi ~ N(0, cr^) by bi ''^ GLG(0,c7, A) with A < 0. The parameter 




-3-2-10 1 2 

Log(baseline period seizures) 



Figure 3: Density of the log(counting) of seizures in the baseline period for the modified data (right). 

estimates, which were obtained by using the procedure NLMIXED in SAS, are described in Table 2. 
Even though the inferential conclusions are the same for both models, indicating a significant effect for 
the treatment group in the sense of decreasing the seizure mean, the random intercept Poisson-GLG 
model seems to fit better the data under the AIC criterion. Furthermore, the estimate A = —1.1852, 
that is significant at 5%, confirms the evidences of Figure 2 on a skew form to right for the random 
intercept distribution. 



13 



Table 2 



Parameter estimates with the respective approximate standard errors for the Poisson-normal 
and Poisson-GLG random intercept models fitted to modified epilepsy data. 







PcnHscni-uoruiMl 






Poissou-GLG 




Parameter 


Estimate 


Sd. error 


z- value 


Estimate 


Sd. error 


z- value 


a 


0.9379 


0.1132 


8.29 


0.5062 


0.2049 


2.47 


P 


0.4841 


0.0418 


11.59 


0.4812 


0.0416 


11.56 


5 


-0.4741 


0.0610 


-7.78 


-0.4685 


0.0606 


-7.72 


a 


0.8426 


0.0802 




0.6130 


0.1418 




A 








-1.1852 


0.6024 


-1.97 






AIC = 2214.62 






AIC = 2155.04 





6.2 C. dubia data 

As a second illustration wc will consider the data set described in Lange et al. (1994) (see also See 
and Bailer, 1988), which was obtained from a reproductive aquatic toxicology experiment and whose 
aim was to study the effect of the herbicide Nitrofen on the asexual reproduction of the freshwater 
invertebrate Ceriodaphnia dubia (C. dubia). The data represent the offspring born counting in three 
broods to each of 10 C. dubia in each of 5 concentration groups of Nitrofen: 0, 80, 160, 235 and 310 
fig /I. The eggs of these species were developed and hatched within 48 hours and were released from 
the brood pouch according to the molting cycle of adult female (24 to 48 hours). The mean profiles 
of the offspring born counting under each concentration across the time are presented in Figure 3. 

Various Poisson regression models with log link were applied by Lange et al. (1994) in order to 
compare the Nitrofen level potencies and the toxin effects on the individual organisms. A Poisson re- 
gression model with a quadratic effect for the concentration was considered as the best-fitting model, 
but the possibility of within-animal correlation (one has three broods for each animal) was not con- 
sidered. The normal probability plot for the deviance residual with generated envelope in Figure 4, 
indicates that the Poisson model is not suitable to fit this data set, with indication of over dispersion. 
Thus, the multivariate negative binomial (MNB) model appears as an option to explain the offspring 
counting and based on the behavior of Figure 3 we suggest the following model to fit the C. dubia 
data: 
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(i) S?MNB(^„(/.) with 

(ii) Hij = exp(/3o + PiCij + ^2Dayj^- + ^aCy x Dayj^), 

where yij is the offspring counting of the ith adult female in the jth brood, for i = 1, ... ,50 and j = 
1, 2, 3, with Cij denoting the concentration for which the jth brood of the ith animal was submitted. 
One has the following settings: Cij = (i = 1,...,10), Cij = 80 (i = 11,..., 20), Q^- = 160 
(i = 21, . . . , 30), Cij = 235 (z = 31, . . . , 40) and Cij = 310 (i = 41, . . . , 50), for j = 1, 2, 3, and Day^j- 
denotes the day in which the eggs were hatched for the jih. brood of the ith animal and it assumes 
the following values: Daya = —2, Dayj2 = and Day^a = 2, for i = 1, . . . , 50. 




Figure 4: Mean profiles of the offspring born counting of C. dubia under each concentration across the 
time. 
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Quantile of theN(0,1) 



Figure 5: Normal probability plot with generated envelope for the Poisson quadratic regression model 
fitted to C. dubia data. 



Table 3 

Parameter estimates with the respective approximate standard errors for the negative binomial 



(XB) 


and mulli\ 


'Mriiilc^ ucgMlive liincnuiid (MXB) 


models lill 


cd lo C. dvbta di\\i 


. 






NB model 






MNB model 




Parameter 


Estimate 


Sd. Error 


z- value 


Estimate 


Sd. Error 


z- value 




2.5269 


0.0704 


35.8950 


2.5333 


0.0925 


27.3619 


Pi 


-0.0040 


0.0004 


-9.8020 


-0.0040 


0.0005 


-7.5045 




0.3329 


0.0432 


7.7050 


0.2916 


0.0297 


9.8088 




-0.0015 


0.0003 


-6.2020 


-0.0013 


0.0002 


-6.8379 




7.4116 


2.1900 




11.5603 


4.2267 




AlC 




830.3 






829 




Deviance 




225.76 (146 d.f.) 






222.91 (146 d.f.) 





The parameter estimates (standard errors) of the NB and MNB models, given in Table 3, are similar, 

confirming the tendencies observed in Figure 3. However, the AIC and deviance values suggest that 
the MNB model fitts better the data. This can also be observed by comparing the normal probability 
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Figure 6: Normal probability plot with generated envelope for the univariate negative binomial model 
(left) and multivariate negative binomial model model (right) fitted to C. dubia data. 

plots for the deviance component and Pearson residuals in Figure 5. Lange et al. affirms that the 
interclass correlations in the C. dubia data is small. For (f) = 11.5603, we identify moderate to weak 
intraclass correlations with the increase of the concentration levels. We also obtain A = 0.293(0.053), 
a since that 4> = When A assumes a small value the MNB model has the feature of fits data sets 
positively correlated with a considerable number of zeros. The data set C. dubia contains many zeros in 
the highest concentration level, however, some of the intraclass correlations are nearly zero. This fact 
can explain the lack of fit observed in Figure 5. The deviance component residual has been suggested 
in this paper was not used because cP {yij , jlij , (j)) assumed any negative values. We are investigating 
this fact. 
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7 Concluding remarks 



In this paper we propose the generahzed log-gamma distribution to give flexibihty for the random 
intercept distribution in Poisson mixed models. The advantage of this distribution is the skew forms 
to the right and to the left including the normal distribution as a particular case. Prom the random 
intercept Poisson-GLG model the multivariate negative binomial model was derived for a particular 
parameter setting and the score functions, Pisher information matrix as well as an iterative process 
were derived. Residual analysis were also proposed. In addition, we present two motivating examples 
emphasizing the specials features of each model. Particularly, for the epilepsy application, we conclude 
that the random intercept Poisson-GLG model seems to be more appropriate to fit the modified data 
with indication of skew form for the intercept distribution as well as presence of outlying observation. 
In the C. dubia application the multivariate negative binomial model seems to fit better the data 
than the univariate negative binomial model, since it incorporates the intraclass correlation that is in 
general positive. Thus, we believe that the models proposed in this work enlarge the options in the 
class of generalized linear mixed models particularly to fit count data with indication of overdispersion 
and nonnormal distribution for the random effects. Extensions for other responses, such binomial and 
gamma, as well as for two or more random effects are in progress. 
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